/**********************************************************************
Compile : gcc -O3 -o generate3 generate3.c
Usage   : ./generate3 N ratio1 ratio2 ratio3 ratio4 M1 M2 seed
	  ./generate3 11 0.3 0.5 0.3 0.3 10000 10 $RANDOM
N       : number of the nodes
ratio1  : number of interactions (>1.0 per node; <1.0 percentage)
ratio2  : inhibition / activation
ratio3  : number of self-degradations
ratio4  : number of activated nodes in the initial state
M1      : the number of the attempts
M2      : the number of the attempts for the same network
seed    : random seed
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

struct netstat {
	int nr; /* nr: number of red edges (inhibition) */
	int ng; /* ng: number of green edges (activations) */
	int ns; /* ns: number of self-degradations */
};

/* state[N] -> state2[N] */
void dynamics (int N, int state[N], int network[N][N], int state2[N])
{
	int j, k, flag;

	for (j = 0; j < N; j++)
	{
		flag = 0;
		for (k = 0; k < N; k++)
		{
			if (k == j)
				continue;
			if (state[k] == 1 && network[j][k] == -1)
			{
				state2[j] = 0;
				flag = 1;
				break;
			}
		}

		if (flag == 1)
			continue;

		for (k = 0; k < N; k++)
		{
			if (state[k] == 1 && network[j][k] == 1)
			{
				state2[j] = 1;
				flag = 1;
				break;
			}
		}

		if (flag == 1)
			continue;

		if (network[j][j] == -1)
			state2[j] = 0;
		else
			state2[j] = state[j];
	}
}

/* Return value : positive fixpoint; negative limited cycle */
int generate_dynamics (int N, int network[N][N], int init[N])
{
	int i, j, n, flag, p;
	struct trajectory {
		int max_L;
		int *state;
	} trajectory;

	trajectory.max_L = 50;
	trajectory.state = calloc (trajectory.max_L, sizeof (int) * N);

	memcpy (trajectory.state, init, sizeof (int) * N);
	n = 1;
	while (1)
	{
		if (n == trajectory.max_L)
		{
			trajectory.max_L += 50;
			trajectory.state = realloc (trajectory.state, trajectory.max_L * sizeof (int) * N);
		}
		dynamics (N, trajectory.state + (n - 1) * N, network, trajectory.state + n * N);
		n ++;

		for (i = n - 1; i > 0; i--)
		{
			flag = 0;
			for (j = 0; j < N; j++)
			{
				if (trajectory.state[(n - 1) * N + j] != trajectory.state[(i - 1) * N + j])
				{
					flag = 1;
					break;
				}
			}
			if (flag == 0)
			{
				break;
			}
		}
		if (flag == 0)
		{
			break;
		}
	}

	if (i == n - 1)
	{
		printf ("Fixpoint : %3d\n", n);
		p = n;
	}
	else
	{
		printf ("Limited cycle : %3d\n", n - i);
		p = i - n;
	}

	/*
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < N; j++)
		{
			printf ("%d ", trajectory.state[i * N + j]);
		}
		printf ("\n");
	}
	*/

	free (trajectory.state);

	return p;
}

/* Return value 1/0: good/bad network */
int check_network (int N, int network[N][N])
{
	int i, j;
	int flag, flag1, flag2;

	flag = 1;
	for (i = 0; i < N; i++)
	{
		flag1 = flag2 = 0;
		for (j = 0; j < N; j++)
		{
			if (network[j][i] == 1)
				flag1 = 1;
			else if (network[j][i] == -1)
				flag2 = 1;
			if (flag1 && flag2)
				break;
		}

		/* make sure every node is regulated by some other nodes */
		if (flag1 == 0 || flag2 == 0)
		{
			flag = 0;
			break;
		}
	}

	for (i = 0; i < N; i++)
	{
		flag1 = 0;
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;
			if (network[i][j] != 0)
			{
				flag1 = 1;
				break;
			}
		}

		/* make sure every node regulated some other nodes */
		if (flag1 == 0)
		{
			flag = 0;
			break;
		}
	}

	return flag;
}

void stat_network (int N, int network[N][N], struct netstat *netstat)
{
	int i, j;

	netstat->nr = netstat->ng = netstat->ns = 0;
	for (i = 0; i < N; i++)
	{
		if (network[i][i] == -1)
			netstat->ns ++;
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;
			if (network[i][j] == -1)
				netstat->nr ++;
			else if (network[i][j] == 1)
				netstat->ng ++;
		}
	}
}

void generate_network (int N, int network[N][N], float ratio1, float ratio2, float ratio3)
{
	int i, j;

	for (i = 0; i < N; i++)
	{
		if (drand48 () < ratio3)
			network[i][i] = -1;
		else
			network[i][i] = 0;

		for (j = 0; j < N; j++)
		{
			if (i == j) continue;
			if (drand48 () < ratio1)
			{
				if (drand48 () < ratio2)
					network[i][j] = -1;
				else
					network[i][j] = 1;
			}
			else
				network[i][j] = 0;
		}
	}
}

void generate_state (int N, int state[N], float ratio)
{
	int i;

	for (i = 0; i < N; i++)
	{
		if (drand48 () < ratio)
			state[i] = 1;
		else
			state[i] = 0;
	}
}

void print_network (int N, int network[N][N])
{
	int i, j;

	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (network[i][j] == 0)
				printf ("n ");
			else if (network[i][j] == 1)
				printf ("g ");
			else if (network[i][j] == -1)
				printf ("r ");
		}
		printf ("\n");
	}
}

int main (int argc, char *argv[])
{
	long int seed;
	int N, M1, M2;
	float ratio1, ratio2, ratio3, ratio4;

	N = atoi (argv[1]);
	ratio1 = atof (argv[2]);
	if (ratio1 > 1.0)
		ratio1 = ratio1 / (float) N;
	ratio2 = atof (argv[3]);
	ratio2 = ratio2 / (ratio2 + 1.0);
	ratio3 = atof (argv[4]);
	if (ratio3 > 1.0)
		ratio3 = ratio3 / (float) N;
	ratio4 = atof (argv[5]);
	if (ratio4 > 1.0)
		ratio4 = ratio4 / (float) N;
	M1 = atoi (argv[6]);
	M2 = atoi (argv[7]);
	seed = atoi (argv[8]);
	srand48 (seed);

	int i, j, flag, L;
	int network[N][N], init[N];
	struct netstat netstat;

	for (i = 0; i < M1; i++)
	{
		generate_network (N, network, ratio1, ratio2, ratio3);
		stat_network (N, network, &netstat);
		flag = check_network (N, network);
		if (flag == 1)
			printf ("c regular network %3d %3d %3d\n", netstat.nr, netstat.ng, netstat.ns);
		else
			printf ("c random network %3d %3d %3d\n", netstat.nr, netstat.ng, netstat.ns);
		//print_network (N, network);
		for (j = 0; j < M2; j++)
		{
			generate_state (N, init, ratio4);
			L = generate_dynamics (N, network, init);
		}
	}
}
